


load(['Mats_T2C.mat'])
Cons_T2C.NETMISS = MISS_T2C.NETMISS;
Cons_T2C.A_MISS = MISS_T2C.A_MISS;
Cons_T2C.y1_educ_MISS = MISS_T2C.y1_educ_MISS;
Cons_T2C.y1_gender_MISS = MISS_T2C.y1_gender_MISS;


y1_raw = struct('y1_educ_m1',NetOut_T2C.y1_educ,'y1_educ_m2',NetOut_T2C.y1_educ,'y1_educ_m3',NetOut_T2C.y1_educ,'y1_educ_m4',NetOut_T2C.y1_educ, ...
    'y1_gender_m1',NetOut_T2C.y1_gender,'y1_gender_m2',NetOut_T2C.y1_gender,'y1_gender_m3',NetOut_T2C.y1_gender,'y1_gender_m4',NetOut_T2C.y1_gender);

y1_imp = struct('y1_educ_m1',NetOut_T2C.y1_educ,'y1_educ_m2',NetOut_T2C.y1_educ,'y1_educ_m3',NetOut_T2C.y1_educ,'y1_educ_m4',NetOut_T2C.y1_educ, ...
    'y1_gender_m1',NetOut_T2C.y1_gender,'y1_gender_m2',NetOut_T2C.y1_gender,'y1_gender_m3',NetOut_T2C.y1_gender,'y1_gender_m4',NetOut_T2C.y1_gender);

G_ij_zero = NetOut_T2C.G_ij;

% To get Var_G and DELTA
M1_G = sum(G_ij_zero.*(ones(size(Cons_T2C.NETMISS,1),1) - Cons_T2C.NETMISS))/sum(ones(size(Cons_T2C.NETMISS,1),1) - Cons_T2C.NETMISS);
M2_G = sum((G_ij_zero.^2).*(ones(size(Cons_T2C.NETMISS,1),1) - Cons_T2C.NETMISS))/sum(ones(size(Cons_T2C.NETMISS,1),1) - Cons_T2C.NETMISS);
V_G = M2_G - M1_G^2;

DELTA = 0;    


for c=1:chains_NoA
    
    c

    Chain_Out_NoA = struct('BGD_hat_b',[],'V_C_b',[],'l_M_hat_b',[],'beta_M_hat_b',[],'sigma2_M_hat_b',[], ...
        'alpha_educ_m1_b',[],'alpha_educ_m3_b',[], ...
        'alpha_gender_m1_b',[],'alpha_gender_m3_b',[], ...
        'y1_educ_m1_b',[],'y1_educ_m3_b',[],...
        'y1_gender_m1_b',[],'y1_gender_m3_b',[], ...
        'Sigma2_educ_m1_b',[],'Sigma2_educ_m3_b',[], ...
        'Sigma2_gender_m1_b',[],'Sigma2_gender_m3_b',[]);    

    G_ij_imp = NetOut_T2C.G_ij + DELTA * Cons_T2C.NETMISS;
    NetEst_output_imp.A_hat = zeros(Cons_T2C.n,1);
        
    %%%%% NETWORK ESTIMATION / IMPUTATION %%%%%
    for w=1:PE_sim_reps
        
        
        
        %%% (1) Estimation Step
        % Before anything, delete DELTA from G_ij_imp
        G_ij_est = G_ij_imp - DELTA * Cons_T2C.NETMISS;
        
        % Reset G_ij's
        G_ji_est = Cons_T2C.Flip_ij*G_ij_est;
        G_ij_doti = Cons_T2C.Trans_doti*G_ij_est;
        %G_ji_doti = Cons_T2C.Trans_doti*Cons_T2C.Flip_ij*G_ij_est;

        % Model 1 (includes X only)
        if Model==1
            % Estimate
            RHS_imp = [G_ji_est, Cons_T2C.X_j(:,1:6), bsxfun(@times,Cons_T2C.X_i(:,1:6),Cons_T2C.X_j(:,1:6))];
            RHS_doti_imp = Cons_T2C.Trans_doti*RHS_imp;
            Z_NoA = [RHS_doti_imp(:,2:size(RHS_doti_imp,2)), Cons_T2C.X_j_mean_Xk_less_ij(:,1:6), Cons_T2C.X_i_mean_Xk_less_ij(:,1:6)];
            BGD_hat_NoA_imp = inv(RHS_doti_imp'*Z_NoA*Z_NoA'*RHS_doti_imp)'*RHS_doti_imp'*Z_NoA*Z_NoA'*G_ij_doti;
            BGD_hat_NoA_imp(1)
            
            % Reset Estimates
            C_hat_raw_imp = -G_ij_imp + RHS_imp*BGD_hat_NoA_imp;
            lambda_hat_imp = inv(Cons_T2C.I_i'*Cons_T2C.I_i)*Cons_T2C.I_i'*C_hat_raw_imp;
            C_ij_hat_imp = C_hat_raw_imp - Cons_T2C.I_i*lambda_hat_imp;
            %C_ji_hat_imp = Cons_T2C.Flip_ij*C_ij_hat_imp;
            
            l_M_hat_imp = log(Cons_T2C.I_i'*(exp(G_ij_imp).*exp(C_ij_hat_imp)));

            X_M = [ones(size(Cons_T2C.X_raw,1),1), Cons_T2C.X_raw];
            beta_M_hat_imp = inv(X_M'*X_M)*X_M'*l_M_hat_imp;
            sigma2_M_hat_imp = mean((l_M_hat_imp - X_M*beta_M_hat_imp).^2);
            
            NetEst_output_imp.BGD_hat = [BGD_hat_NoA_imp(1:7); zeros(7,1); 0; BGD_hat_NoA_imp(8:13); zeros(34,1)];
            %NetEst_output_imp.C_ij_hat = C_ij_hat_imp;
            NetEst_output_imp.l_M_hat = l_M_hat_imp;
            NetEst_output_imp.beta_M_hat = [beta_M_hat_imp; zeros(7,1)];
            NetEst_output_imp.sigma2_M_hat = sigma2_M_hat_imp;
            NetEst_output_imp.V_C = [C_ij_hat_imp, Cons_T2C.Flip_ij*C_ij_hat_imp]'*[C_ij_hat_imp, Cons_T2C.Flip_ij*C_ij_hat_imp]/size(C_ij_hat_imp,1);               
        end
        
        % Model 2 (includes X, P, and interactions)
        if Model==2
            RHS_imp = [G_ji_est, Cons_T2C.X_j, bsxfun(@times,Cons_T2C.X_i,Cons_T2C.X_j)];
            RHS_doti_imp = Cons_T2C.Trans_doti*RHS_imp;
            Z_NoA = [RHS_doti_imp(:,2:size(RHS_doti_imp,2)), Cons_T2C.X_j_mean_Xk_less_ij, Cons_T2C.X_i_mean_Xk_less_ij];
            BGD_hat_NoA_imp = inv(RHS_doti_imp'*Z_NoA*Z_NoA'*RHS_doti_imp)'*RHS_doti_imp'*Z_NoA*Z_NoA'*G_ij_doti;
            BGD_hat_NoA_imp(1)
            
            % Reset Estimates
            C_hat_raw_imp = -G_ij_imp + RHS_imp*BGD_hat_NoA_imp;
            lambda_hat_imp = inv(Cons_T2C.I_i'*Cons_T2C.I_i)*Cons_T2C.I_i'*C_hat_raw_imp;
            C_ij_hat_imp = C_hat_raw_imp - Cons_T2C.I_i*lambda_hat_imp;
            %C_ji_hat_imp = Cons_T2C.Flip_ij*C_ij_hat_imp;
            
            l_M_hat_imp = log(Cons_T2C.I_i'*(exp(G_ij_imp).*exp(C_ij_hat_imp)));

            X_M = [ones(size(Cons_T2C.X_raw,1),1), Cons_T2C.X_raw, Cons_T2C.P, bsxfun(@times,Cons_T2C.P,Cons_T2C.X_raw)];
            beta_M_hat_imp = inv(X_M'*X_M)*X_M'*l_M_hat_imp;
            sigma2_M_hat_imp = mean((l_M_hat_imp - X_M*beta_M_hat_imp).^2);   



            NetEst_output_imp.BGD_hat = [BGD_hat_NoA_imp(1:14); 0; BGD_hat_NoA_imp(15:27); zeros(27,1)];
            %NetEst_output_imp.C_ij_hat = C_ij_hat_imp;
            NetEst_output_imp.l_M_hat = l_M_hat_imp;
            NetEst_output_imp.beta_M_hat = beta_M_hat_imp;
            NetEst_output_imp.sigma2_M_hat = sigma2_M_hat_imp;
            NetEst_output_imp.V_C = [C_ij_hat_imp, Cons_T2C.Flip_ij*C_ij_hat_imp]'*[C_ij_hat_imp, Cons_T2C.Flip_ij*C_ij_hat_imp]/size(C_ij_hat_imp,1);            
        end


        %%% (2) Imputation Step
        [G_ij_imp_new, A_imp_new] = Fn_ImputeNetwork_Mar2022(G_ij_imp, NetEst_output_imp, Cons_T2C, sim_tol, DELTA);

        % Reset Estimates                                
        G_ij_imp = G_ij_imp_new;

             
    end
    
    V_C_imp = [C_ij_hat_imp, Cons_T2C.Flip_ij*C_ij_hat_imp]'*[C_ij_hat_imp, Cons_T2C.Flip_ij*C_ij_hat_imp]/size(C_ij_hat_imp,1);
    

    %%%%% PEER EFFECTS ESTIMATION / IMPUTATION %%%%%
    if Net_OUT==0
        W_ij = exp(G_ij_imp) + exp(Cons_T2C.Flip_ij*G_ij_imp);
    else
        W_ij = exp(G_ij_imp);
    end    
    
    W_ij = Cons_T2C.W_trans*W_ij;
    W_ij_rep = zeros(Cons_T2C.n,Cons_T2C.n);
    k = 0;
    for s=1:Cons_T2C.schools
        W_ij_rep(k+1:k+Cons_T2C.sch_size(s),k+1:k+Cons_T2C.sch_size(s)) = reshape(W_ij(k+1:k+Cons_T2C.sch_size(s)^2,1),Cons_T2C.sch_size(s),Cons_T2C.sch_size(s));
        k = k + Cons_T2C.sch_size(s);
    end
    W_ij_imp = bsxfun(@rdivide,W_ij_rep,sum(W_ij_rep,2));   

    A_hat = zeros(Cons_T2C.n,1);
    

    for w=1:PE_sim_reps
        w
        %% Estimation(A_hat, W_ij, y1, Cons)
        [PE_out_imp] = Fn_EstimatePE_Mar2022(A_hat, W_ij_imp, y1_imp, Cons_T2C); 
        
        %% Imputation Step
        [y1_imp] = Fn_ImputePE_Mar2022(PE_out_imp, A_hat, W_ij_imp, Cons_T2C, y1_raw);

    end    
    
    %% Save estimates and imputations
    Chain_Out_NoA.BGD_hat_b = [Chain_Out_NoA.BGD_hat_b, BGD_hat_NoA_imp];
    %Chain_Out_NoA.G_ij_b = [Chain_Out_NoA.G_ij_b, G_ij_imp];
    %Chain_Out_NoA.C_ij_hat_b = [Chain_Out_NoA.C_ij_hat_b, C_ij_hat_imp];
    Chain_Out_NoA.V_C_b = [Chain_Out_NoA.V_C_b, V_C_imp(:,1)];
    Chain_Out_NoA.l_M_hat_b = [Chain_Out_NoA.l_M_hat_b, l_M_hat_imp];
    Chain_Out_NoA.beta_M_hat_b = [Chain_Out_NoA.beta_M_hat_b, beta_M_hat_imp];
    Chain_Out_NoA.sigma2_M_hat_b = [Chain_Out_NoA.sigma2_M_hat_b, sigma2_M_hat_imp];
    
    
    %% For estimates, only save models 1 and 3 (2 and 4 involve A and are thus undefined
    % Estimates
    Chain_Out_NoA.alpha_educ_m1_b = [Chain_Out_NoA.alpha_educ_m1_b, PE_out_imp.alpha_educ_m1];
    Chain_Out_NoA.alpha_educ_m3_b = [Chain_Out_NoA.alpha_educ_m3_b, PE_out_imp.alpha_educ_m3];
    Chain_Out_NoA.alpha_gender_m1_b = [Chain_Out_NoA.alpha_gender_m1_b, PE_out_imp.alpha_gender_m1];
    Chain_Out_NoA.alpha_gender_m3_b = [Chain_Out_NoA.alpha_gender_m3_b, PE_out_imp.alpha_gender_m3];

    % imputed Values
    Chain_Out_NoA.y1_educ_m1_b = [Chain_Out_NoA.y1_educ_m1_b, y1_imp.y1_educ_m1];
    Chain_Out_NoA.y1_educ_m3_b = [Chain_Out_NoA.y1_educ_m3_b, y1_imp.y1_educ_m3];
    Chain_Out_NoA.y1_gender_m1_b = [Chain_Out_NoA.y1_gender_m1_b, y1_imp.y1_gender_m1];
    Chain_Out_NoA.y1_gender_m3_b = [Chain_Out_NoA.y1_gender_m3_b, y1_imp.y1_gender_m3];
    
    % Implied error variances    
    Chain_Out_NoA.Sigma2_educ_m1_b = [Chain_Out_NoA.Sigma2_educ_m1_b, PE_out_imp.Sigma2_educ_m1];
    Chain_Out_NoA.Sigma2_educ_m3_b = [Chain_Out_NoA.Sigma2_educ_m3_b, PE_out_imp.Sigma2_educ_m3];
    Chain_Out_NoA.Sigma2_gender_m1_b = [Chain_Out_NoA.Sigma2_gender_m1_b, PE_out_imp.Sigma2_gender_m1];
    Chain_Out_NoA.Sigma2_gender_m3_b = [Chain_Out_NoA.Sigma2_gender_m3_b, PE_out_imp.Sigma2_gender_m3];



end

if Net_OUT==0
    save(['Estimates_NoA_M',num2str(Model),'_OR.mat'],'Chain_Out_NoA')
else
    save(['Estimates_NoA_M',num2str(Model),'_OUT.mat'],'Chain_Out_NoA')
end
    
    
